###########################
### PANDEMICS AND POLITICAL DEVELOPMENT
### THE ELECTORAL LEGACY OF THE BLACK DEATH IN GERMANY
### WORLD POLITICS, VOL. 73, NO. 3 (JULY 2021)
### REPLICATION FILE OF THE EMPIRICAL ANALYSIS
### DANIEL W. GINGERICH & JAN P. VOGLER
###########################

###########################
### REPLICATION FILES PART 2
### WEIMAR GERMANY
###########################

###########################
### INSTALL REQUIRED PACKAGES
###########################

install.packages("foreign")
install.packages("geosphere")
install.packages("AER")
install.packages("multiwayvcov")
install.packages("lmtest")
install.packages("ggplot2")
install.packages("egg")



###########################
### DATA MANAGEMENT
###########################

###########################
### DATA ON ELECTORAL OUTCOMES IN WEIMAR GERMANY
### BY FALTER & HAENISCH (1990)
###########################

library(foreign)

data_outcomes=read.dta("ZA8013_Wahldaten.dta")

# CHANGE "-9" (NA) ENTRIES TO "NA"

data_outcomes$n309gs[which(data_outcomes$n309gs==-9)]=NA # VALID VOTES, 1930 ELECTION
data_outcomes$n309nsda[which(data_outcomes$n309nsda==-9)]=NA # VOTES FOR THE NSDAP, 1930 ELECTION
data_outcomes$n309dnvp[which(data_outcomes$n309dnvp==-9)]=NA # VOTES FOR THE DNVP, 1930 ELECTION

data_outcomes$n327gs[which(data_outcomes$n327gs==-9)]=NA # VALID VOTES, JULY 1932 ELECTION
data_outcomes$n327nsda[which(data_outcomes$n327nsda==-9)]=NA # VOTES FOR THE NSDAP, JULY 1932 ELECTION
data_outcomes$n327dnvp[which(data_outcomes$n309dnvp==-9)]=NA # VOTES FOR THE DNVP, JULY 1932 ELECTION

# CREATE OUTCOME VARIABLES

data_outcomes$nsdap1930=data_outcomes$n309nsda/data_outcomes$n309gs # NSDAP VOTE SHARE (1930)
data_outcomes$nsdap1932j=data_outcomes$n327nsda/data_outcomes$n327gs # NSDAP VOTE SHARE (JUL. 1932)
data_outcomes$rightwing1930=(data_outcomes$n309nsda+data_outcomes$n309dnvp)/data_outcomes$n309gs # NSDAP AND DNVP COMBINED VOTE SHARE (1930)
data_outcomes$rightwing1932j=(data_outcomes$n327nsda+data_outcomes$n327dnvp)/data_outcomes$n327gs # NSDAP AND DNVP COMBINED VOTE SHARE (JUL. 1932)



###########################
### DATA ON CONTROL VARIABLES
###########################

data_control_var=read.csv("data_wei_control_var.csv", stringsAsFactors=F)



###########################
### DATA ON TOWN/CITY/COUNTY GEOLOCATION
### BY SELB & MUNZERT (2018)
###########################

data_weimar_geolocation=read.delim("geopoints_df_cleaned.csv", header=T, sep=";", quote="", fill=FALSE)

data_weimar_geolocation=data_weimar_geolocation[,c("krnr","lfnr","lat","long")]



###########################
### DATA ON BLACK DEATH OUTBREAKS
### BY JEDWAB ET AL. (2019)
###########################

data_bd_outbreaks=read.csv("data_wei_bd_outbreaks.csv", stringsAsFactors=F)



###########################
### MERGE DATASETS FOR ANALYSIS
###########################

data_weimar_fin=merge(data_outcomes, data_weimar_geolocation, by=c("krnr","lfnr"), all.x=T)

data_weimar_fin=merge(data_weimar_fin, data_control_var, by=c("krnr","lfnr"), all.x=T)



###########################
### CALCULATING BDEI SCORES
###########################

###########################
### COMPUTE THE GLOBAL MAXIMUM DISTANCE
###########################

library(geosphere)
max_store=rep(NA,length(data_weimar_fin$lfnr))
for (j in 1:length(data_weimar_fin$lfnr)){
  store=rep(NA,length(data_bd_outbreaks$int))
  for (i in 1:length(data_bd_outbreaks$int)){
    store[i]=distm(c(data_weimar_fin$long[j], data_weimar_fin$lat[j]), c(data_bd_outbreaks$longitude[i], data_bd_outbreaks$latitude[i]), fun = distHaversine)
    max_store[j]=max(store)
  }
}
max_dist=max(max_store)



###########################
### COMPUTE BLACK DEATH EXPOSURE INTENSITY SCORE (V1-5)
###########################

for (j in 1:length(data_weimar_fin$lfnr)){
  store=rep(NA,length(data_bd_outbreaks$int))
  for (i in 1:length(data_bd_outbreaks$int)){
    store[i]=distm(c(data_weimar_fin$long[j], data_weimar_fin$lat[j]), c(data_bd_outbreaks$longitude[i], data_bd_outbreaks$latitude[i]), fun = distHaversine)
    data_weimar_fin$BDEI_score[j]=sum(data_bd_outbreaks$int*(1-store/max_dist)^3)
    data_weimar_fin$BDEI_score2[j]=sum(data_bd_outbreaks$int*(1-store/max_dist)^6)
    data_weimar_fin$BDEI_score3[j]=sum(data_bd_outbreaks$int*(1-store/max_dist)^9)
    data_weimar_fin$BDEI_score4[j]=sum(data_bd_outbreaks$int*(1-store/max_dist)^12)
    data_weimar_fin$BDEI_score5[j]=sum(data_bd_outbreaks$int*(1-store/max_dist)^15)
  }
}



###########################
### STANDARDIZE BLACK DEATH EXPOSURE INTENSITY SCORE (V1-5)
###########################

data_weimar_fin$BDEI_score=scale(data_weimar_fin$BDEI_score)[,1]
data_weimar_fin$BDEI_score2=scale(data_weimar_fin$BDEI_score2)[,1]
data_weimar_fin$BDEI_score3=scale(data_weimar_fin$BDEI_score3)[,1]
data_weimar_fin$BDEI_score4=scale(data_weimar_fin$BDEI_score4)[,1]
data_weimar_fin$BDEI_score5=scale(data_weimar_fin$BDEI_score5)[,1]



###########################
### CREATE SUBSET FOR 1930 ANALYSIS
###########################

data_weimar_fin_analysis=data_weimar_fin
data_weimar_fin_analysis=data_weimar_fin_analysis[which(data_weimar_fin$agglvl=="GEMEINDEN AB 2000 E." | data_weimar_fin$agglvl=="stadtkreise"),]



###########################
### CREATE SUBSET FOR JUL. 1932 ANALYSIS
###########################

data_weimar_fin_1932_analysis=data_weimar_fin
data_weimar_fin_1932_analysis=data_weimar_fin_1932_analysis[-which(is.na(data_weimar_fin_1932_analysis$nsdap1932j)),]



###########################
### MAIN EMPIRICAL ANALYSIS IN THE ARTICLE
###########################

###########################
### NSDAP VOTE SHARE (1930)
### TABLE 4 IN THE ARTICLE
###########################

vote1930_lm1=lm(nsdap1930 ~ BDEI_score, data=data_weimar_fin_analysis)
summary(vote1930_lm1)

vote1930_lm2=lm(nsdap1930 ~ BDEI_score2, data=data_weimar_fin_analysis)
summary(vote1930_lm2)

vote1930_lm3=lm(nsdap1930 ~ BDEI_score3, data=data_weimar_fin_analysis)
summary(vote1930_lm3)

vote1930_lm4=lm(nsdap1930 ~ BDEI_score4, data=data_weimar_fin_analysis)
summary(vote1930_lm4)

vote1930_lm5=lm(nsdap1930 ~ BDEI_score5, data=data_weimar_fin_analysis)
summary(vote1930_lm5)

vote1930_control_lm1=lm(nsdap1930 ~ BDEI_score + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_analysis)
summary(vote1930_control_lm1)

vote1930_control_lm2=lm(nsdap1930 ~ BDEI_score2 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_analysis)
summary(vote1930_control_lm2)

vote1930_control_lm3=lm(nsdap1930 ~ BDEI_score3 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_analysis)
summary(vote1930_control_lm3)

vote1930_control_lm4=lm(nsdap1930 ~ BDEI_score4 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_analysis)
summary(vote1930_control_lm4)

vote1930_control_lm5=lm(nsdap1930 ~ BDEI_score5 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_analysis)
summary(vote1930_control_lm5)

### CLUSTERED STANDARD ERRORS

library(multiwayvcov)
library(lmtest)

vote1930_lm1.vcovCL=cluster.vcov(vote1930_lm1, data_weimar_fin_analysis$wkr)
vote1930_lm1.clse=coeftest(vote1930_lm1, vote1930_lm1.vcovCL)
vote1930_lm1.clse

vote1930_lm2.vcovCL=cluster.vcov(vote1930_lm2, data_weimar_fin_analysis$wkr)
vote1930_lm2.clse=coeftest(vote1930_lm2, vote1930_lm2.vcovCL)
vote1930_lm2.clse

vote1930_lm3.vcovCL=cluster.vcov(vote1930_lm3, data_weimar_fin_analysis$wkr)
vote1930_lm3.clse=coeftest(vote1930_lm3, vote1930_lm3.vcovCL)
vote1930_lm3.clse

vote1930_lm4.vcovCL=cluster.vcov(vote1930_lm4, data_weimar_fin_analysis$wkr)
vote1930_lm4.clse=coeftest(vote1930_lm4, vote1930_lm4.vcovCL)
vote1930_lm4.clse

vote1930_lm5.vcovCL=cluster.vcov(vote1930_lm5, data_weimar_fin_analysis$wkr)
vote1930_lm5.clse=coeftest(vote1930_lm5, vote1930_lm5.vcovCL)
vote1930_lm5.clse

vote1930_control_lm1.vcovCL=cluster.vcov(vote1930_control_lm1, data_weimar_fin_analysis$wkr)
vote1930_control_lm1.clse=coeftest(vote1930_control_lm1, vote1930_control_lm1.vcovCL)
vote1930_control_lm1.clse

vote1930_control_lm2.vcovCL=cluster.vcov(vote1930_control_lm2, data_weimar_fin_analysis$wkr)
vote1930_control_lm2.clse=coeftest(vote1930_control_lm2, vote1930_control_lm2.vcovCL)
vote1930_control_lm2.clse

vote1930_control_lm3.vcovCL=cluster.vcov(vote1930_control_lm3, data_weimar_fin_analysis$wkr)
vote1930_control_lm3.clse=coeftest(vote1930_control_lm3, vote1930_control_lm3.vcovCL)
vote1930_control_lm3.clse

vote1930_control_lm4.vcovCL=cluster.vcov(vote1930_control_lm4, data_weimar_fin_analysis$wkr)
vote1930_control_lm4.clse=coeftest(vote1930_control_lm4, vote1930_control_lm4.vcovCL)
vote1930_control_lm4.clse

vote1930_control_lm5.vcovCL=cluster.vcov(vote1930_control_lm5, data_weimar_fin_analysis$wkr)
vote1930_control_lm5.clse=coeftest(vote1930_control_lm5, vote1930_control_lm5.vcovCL)
vote1930_control_lm5.clse



###########################
### NSDAP VOTE SHARE (JUL. 1932)
### TABLE 5 IN THE ARTICLE
###########################

vote1932_lm1=lm(nsdap1932j ~ BDEI_score, data=data_weimar_fin_1932_analysis)
summary(vote1932_lm1)

vote1932_lm2=lm(nsdap1932j ~ BDEI_score2, data=data_weimar_fin_1932_analysis)
summary(vote1932_lm2)

vote1932_lm3=lm(nsdap1932j ~ BDEI_score3, data=data_weimar_fin_1932_analysis)
summary(vote1932_lm3)

vote1932_lm4=lm(nsdap1932j ~ BDEI_score4, data=data_weimar_fin_1932_analysis)
summary(vote1932_lm4)

vote1932_lm5=lm(nsdap1932j ~ BDEI_score5, data=data_weimar_fin_1932_analysis)
summary(vote1932_lm5)

vote1932_control_lm1=lm(nsdap1932j ~ BDEI_score + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_1932_analysis)
summary(vote1932_control_lm1)

vote1932_control_lm2=lm(nsdap1932j ~ BDEI_score2 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_1932_analysis)
summary(vote1932_control_lm2)

vote1932_control_lm3=lm(nsdap1932j ~ BDEI_score3 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_1932_analysis)
summary(vote1932_control_lm3)

vote1932_control_lm4=lm(nsdap1932j ~ BDEI_score4 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_1932_analysis)
summary(vote1932_control_lm4)

vote1932_control_lm5=lm(nsdap1932j ~ BDEI_score5 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_1932_analysis)
summary(vote1932_control_lm5)

### CLUSTERED STANDARD ERRORS

vote1932_lm1.vcovCL=cluster.vcov(vote1932_lm1, data_weimar_fin_1932_analysis$wkr)
vote1932_lm1.clse=coeftest(vote1932_lm1, vote1932_lm1.vcovCL)
vote1932_lm1.clse

vote1932_lm2.vcovCL=cluster.vcov(vote1932_lm2, data_weimar_fin_1932_analysis$wkr)
vote1932_lm2.clse=coeftest(vote1932_lm2, vote1932_lm2.vcovCL)
vote1932_lm2.clse

vote1932_lm3.vcovCL=cluster.vcov(vote1932_lm3, data_weimar_fin_1932_analysis$wkr)
vote1932_lm3.clse=coeftest(vote1932_lm3, vote1932_lm3.vcovCL)
vote1932_lm3.clse

vote1932_lm4.vcovCL=cluster.vcov(vote1932_lm4, data_weimar_fin_1932_analysis$wkr)
vote1932_lm4.clse=coeftest(vote1932_lm4, vote1932_lm4.vcovCL)
vote1932_lm4.clse

vote1932_lm5.vcovCL=cluster.vcov(vote1932_lm5, data_weimar_fin_1932_analysis$wkr)
vote1932_lm5.clse=coeftest(vote1932_lm5, vote1932_lm5.vcovCL)
vote1932_lm5.clse

vote1932_control_lm1.vcovCL=cluster.vcov(vote1932_control_lm1, data_weimar_fin_1932_analysis$wkr)
vote1932_control_lm1.clse=coeftest(vote1932_control_lm1, vote1932_control_lm1.vcovCL)
vote1932_control_lm1.clse

vote1932_control_lm2.vcovCL=cluster.vcov(vote1932_control_lm2, data_weimar_fin_1932_analysis$wkr)
vote1932_control_lm2.clse=coeftest(vote1932_control_lm2, vote1932_control_lm2.vcovCL)
vote1932_control_lm2.clse

vote1932_control_lm3.vcovCL=cluster.vcov(vote1932_control_lm3, data_weimar_fin_1932_analysis$wkr)
vote1932_control_lm3.clse=coeftest(vote1932_control_lm3, vote1932_control_lm3.vcovCL)
vote1932_control_lm3.clse

vote1932_control_lm4.vcovCL=cluster.vcov(vote1932_control_lm4, data_weimar_fin_1932_analysis$wkr)
vote1932_control_lm4.clse=coeftest(vote1932_control_lm4, vote1932_control_lm4.vcovCL)
vote1932_control_lm4.clse

vote1932_control_lm5.vcovCL=cluster.vcov(vote1932_control_lm5, data_weimar_fin_1932_analysis$wkr)
vote1932_control_lm5.clse=coeftest(vote1932_control_lm5, vote1932_control_lm5.vcovCL)
vote1932_control_lm5.clse



###########################
### FIGURES IN THE ARTICLE
###########################

###########################
### PREDICT FUNCTION WITH ADJUSTED VARIANCE-COVARIANCE MATRIX
###########################

### BASED ON: https://stackoverflow.com/questions/3790116/using-clustered-covariance-matrix-in-predict-lm

predict.rob = function(x,clcov,newdata){
  if(missing(newdata)){ newdata = x$model }
  tt = terms(x)
  Terms = delete.response(tt)
  m.mat = model.matrix(Terms,data=newdata)
  m.coef = x$coef
  fit = as.vector(m.mat %*% x$coef)
  se.fit = sqrt(diag(m.mat%*%clcov%*%t(m.mat)))
  return(list(fit=fit,se.fit=se.fit))}



###########################
### PREDICTED VALUE PLOTS
###########################

###########################
#### CREATE FIGURE: BDEI SCORE V1 AND NSDAP VOTE SHARE (1930)
###########################

###########################
#### NEW DATA FRAME
###########################

nd = data.frame(BDEI_score=seq(min(data_weimar_fin_analysis$BDEI_score, na.rm=T), max(data_weimar_fin_analysis$BDEI_score, na.rm=T),length.out=10))



###########################
#### PREDICT OUTCOME VALUES
###########################

pred.p1 = predict.rob(vote1930_lm1, clcov=vote1930_lm1.vcovCL, newdata=nd)
pred.p1.table = cbind(pred.p1$fit, pred.p1$se.fit)
fit.p1 = pred.p1
fit.p1 = pred.p1$fit
low.p1 = pred.p1$fit - 1.96*pred.p1$se.fit
high.p1 = pred.p1$fit + 1.96*pred.p1$se.fit
cis.p1 = cbind(fit.p1, low.p1, high.p1)



###########################
### CREATE TWO SEPARATE PLOTS
### PLOT 1: PREDICTED VALUES
### PLOT 2: HISTOGRAM OF INDEPENDENT VARIABLE
###########################

library(ggplot2)

pred.p1.gg=as.data.frame(pred.p1)
cis.p1.gg=as.data.frame(cis.p1)

plot1=ggplot(pred.p1.gg, aes(x=nd$BDEI_score, y=fit.p1)) +
  geom_ribbon(aes(ymin=cis.p1.gg$low.p1, ymax=cis.p1.gg$high.p1), alpha=0.1, fill="Blue", color="Blue", linetype = 2) + scale_fill_brewer(palette="Blues") + ylab("NSDAP Vote Share (1930)") + ggtitle("Predicted Values Plot: BDEI Score v1 and NSDAP Vote Share (1930)") + geom_line(size=2) + theme(axis.text=element_text(size=16,face="bold",color="Black"), axis.title=element_text(size=16,face="bold"), title=element_text(size=16,face="bold"), legend.text=element_text(size=16,face="bold",color="Black")) + scale_x_continuous(limits=c(min(nd$BDEI_score), max(nd$BDEI_score))) +
  theme(plot.title = element_text(hjust = 0.5), axis.title.x = element_blank())

plot1

plot2=ggplot(data_weimar_fin_analysis, aes(x=data_weimar_fin_analysis$BDEI_score)) + geom_histogram(aes(BDEI_score), position="identity", linetype=1, fill="gray60", data = data_weimar_fin_analysis, alpha=0.5, bins = 30) + scale_y_continuous(expand = c(0,0)) + scale_x_continuous(limits=c(min(nd$BDEI_score), max(nd$BDEI_score))) + xlab("BDEI Score v1 (Min-Max)") + ylab("Count") + theme(axis.text=element_text(size=16,face="bold",color="Black"), axis.title=element_text(size=16,face="bold"), title=element_text(size=16,face="bold"), legend.text=element_text(size=16,face="bold",color="Black"))

plot2



###########################
### COMBINE PLOT 1 AND PLOT 2 IN A SINGLE FIGURE
### FIGURE 7 IN THE ARTICLE
###########################

library(egg)
egg::ggarrange(plot1, plot2, heights = c(0.75, 0.25))



###########################
#### CREATE FIGURE: BDEI SCORE V1 AND NSDAP VOTE SHARE (JUL. 1932)
###########################

###########################
### NEW DATA FRAME
###########################

nd2 = data.frame(BDEI_score=seq(min(data_weimar_fin_1932_analysis$BDEI_score, na.rm=T), max(data_weimar_fin_1932_analysis$BDEI_score, na.rm=T),length.out=10))



###########################
#### PREDICT OUTCOME VALUES
###########################

pred.p2 = predict.rob(vote1932_lm1, clcov=vote1932_lm1.vcovCL, newdata=nd2)
pred.p2.table = cbind(pred.p2$fit, pred.p2$se.fit)
fit.p2 = pred.p2
fit.p2 = pred.p2$fit
low.p2 = pred.p2$fit - 1.96*pred.p2$se.fit
high.p2 = pred.p2$fit + 1.96*pred.p2$se.fit
cis.p2 = cbind(fit.p2, low.p2, high.p2)



###########################
### CREATE TWO SEPARATE PLOTS
### PLOT 1: PREDICTED VALUES
### PLOT 2: HISTOGRAM OF INDEPENDENT VARIABLE
###########################

pred.p2.gg=as.data.frame(pred.p2)
cis.p2.gg=as.data.frame(cis.p2)

plot1=ggplot(pred.p2.gg, aes(x=nd$BDEI_score, y=fit)) +
  geom_ribbon(aes(ymin=cis.p2.gg$low.p2, ymax=cis.p2.gg$high.p2), alpha=0.1, fill="Blue", color="Blue", linetype = 2) + scale_fill_brewer(palette="Blues") + xlab("BDEI Score v1 (Min-Max)") + ylab("NSDAP Vote Share (1932)") + ggtitle("Predicted Values Plot: BDEI Score v1 and NSDAP Vote Share (Jul. 1932)") + geom_line(size=2) + theme(axis.text=element_text(size=16,face="bold",color="Black"), axis.title=element_text(size=16,face="bold"), title=element_text(size=16,face="bold"), legend.text=element_text(size=16,face="bold",color="Black")) +
  theme(plot.title = element_text(hjust = 0.5), axis.title.x = element_blank()) + scale_x_continuous(limits=c(min(nd$BDEI_score), max(nd$BDEI_score)))

plot1

plot2=ggplot(data_weimar_fin_1932_analysis, aes(x=data_weimar_fin_1932_analysis$BDEI_score)) + geom_histogram(aes(BDEI_score), position="identity", linetype=1, fill="gray60", data = data_weimar_fin_1932_analysis, alpha=0.5, bins = 30) + scale_y_continuous(expand = c(0,0)) + scale_x_continuous(limits=c(min(nd$BDEI_score), max(nd$BDEI_score))) + xlab("BDEI Score v1 (Min-Max)") + ylab("Count") + theme(axis.text=element_text(size=16,face="bold",color="Black"), axis.title=element_text(size=16,face="bold"), title=element_text(size=16,face="bold"), legend.text=element_text(size=16,face="bold",color="Black"))

plot2



###########################
### COMBINE PLOT 1 AND PLOT 2 IN A SINGLE FIGURE
### FIGURE 8 IN THE ARTICLE
###########################

egg::ggarrange(plot1, plot2, heights = c(0.75, 0.25))



###########################
### ADDITIONAL EMPIRICAL ANALYSIS IN THE SUPPLEMENTARY MATERIAL
###########################

###########################
### DESCRIPTIVE SUMMARY STATISTICS
### TABLE A21 IN THE SUPPLEMENTARY MATERIAL
###########################

summary(data_weimar_fin)



###########################
### TOBIT MODELS AS AN ALTERNATIVE SPECIFICATION
###########################

###########################
### NSDAP VOTE SHARE (1930)
### TABLE A22 IN THE SUPPLEMENTARY MATERIAL
###########################

library(AER)

vote1930_tobit1=tobit(nsdap1930 ~ BDEI_score, left=0, right=1, data=data_weimar_fin_analysis)
summary(vote1930_tobit1)

vote1930_tobit2=tobit(nsdap1930 ~ BDEI_score2, left=0, right=1, data=data_weimar_fin_analysis)
summary(vote1930_tobit2)

vote1930_tobit3=tobit(nsdap1930 ~ BDEI_score3, left=0, right=1, data=data_weimar_fin_analysis)
summary(vote1930_tobit3)

vote1930_tobit4=tobit(nsdap1930 ~ BDEI_score4, left=0, right=1, data=data_weimar_fin_analysis)
summary(vote1930_tobit4)

vote1930_tobit5=tobit(nsdap1930 ~ BDEI_score5, left=0, right=1, data=data_weimar_fin_analysis)
summary(vote1930_tobit5)

vote1930_control_tobit1=tobit(nsdap1930 ~ BDEI_score + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=100, data=data_weimar_fin_analysis)
summary(vote1930_control_tobit1)

vote1930_control_tobit2=tobit(nsdap1930 ~ BDEI_score2 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=100, data=data_weimar_fin_analysis)
summary(vote1930_control_tobit2)

vote1930_control_tobit3=tobit(nsdap1930 ~ BDEI_score3 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=100, data=data_weimar_fin_analysis)
summary(vote1930_control_tobit3)

vote1930_control_tobit4=tobit(nsdap1930 ~ BDEI_score4 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=100, data=data_weimar_fin_analysis)
summary(vote1930_control_tobit4)

vote1930_control_tobit5=tobit(nsdap1930 ~ BDEI_score5 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=100, data=data_weimar_fin_analysis)
summary(vote1930_control_tobit5)



###########################
### NSDAP VOTE SHARE (JUL. 1932)
### TABLE A23 IN THE SUPPLEMENTARY MATERIAL
###########################

vote1932_tobit1=tobit(nsdap1932j ~ BDEI_score, left=0, right=1, data=data_weimar_fin_1932_analysis)
summary(vote1932_tobit1)

vote1932_tobit2=tobit(nsdap1932j ~ BDEI_score2, left=0, right=1, data=data_weimar_fin_1932_analysis)
summary(vote1932_tobit2)

vote1932_tobit3=tobit(nsdap1932j ~ BDEI_score3, left=0, right=1, data=data_weimar_fin_1932_analysis)
summary(vote1932_tobit3)

vote1932_tobit4=tobit(nsdap1932j ~ BDEI_score4, left=0, right=1, data=data_weimar_fin_1932_analysis)
summary(vote1932_tobit4)

vote1932_tobit5=tobit(nsdap1932j ~ BDEI_score5, left=0, right=1, data=data_weimar_fin_1932_analysis)
summary(vote1932_tobit5)

vote1932_control_tobit1=tobit(nsdap1932j ~ BDEI_score + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=100, data=data_weimar_fin_1932_analysis)
summary(vote1932_control_tobit1)

vote1932_control_tobit2=tobit(nsdap1932j ~ BDEI_score2 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=100, data=data_weimar_fin_1932_analysis)
summary(vote1932_control_tobit2)

vote1932_control_tobit3=tobit(nsdap1932j ~ BDEI_score3 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=100, data=data_weimar_fin_1932_analysis)
summary(vote1932_control_tobit3)

vote1932_control_tobit4=tobit(nsdap1932j ~ BDEI_score4 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=100, data=data_weimar_fin_1932_analysis)
summary(vote1932_control_tobit4)

vote1932_control_tobit5=tobit(nsdap1932j ~ BDEI_score5 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=100, data=data_weimar_fin_1932_analysis)
summary(vote1932_control_tobit5)



###########################
### NSDAP AND DNVP COMBINED VOTE SHARE (1930)
### TABLE A24 IN THE SUPPLEMENTARY MATERIAL
###########################

rightwing1930_lm1=lm(rightwing1930 ~ BDEI_score, data=data_weimar_fin_analysis)
summary(rightwing1930_lm1)

rightwing1930_lm2=lm(rightwing1930 ~ BDEI_score2, data=data_weimar_fin_analysis)
summary(rightwing1930_lm2)

rightwing1930_lm3=lm(rightwing1930 ~ BDEI_score3, data=data_weimar_fin_analysis)
summary(rightwing1930_lm3)

rightwing1930_lm4=lm(rightwing1930 ~ BDEI_score4, data=data_weimar_fin_analysis)
summary(rightwing1930_lm4)

rightwing1930_lm5=lm(rightwing1930 ~ BDEI_score5, data=data_weimar_fin_analysis)
summary(rightwing1930_lm5)

rightwing1930_control_lm1=lm(rightwing1930 ~ BDEI_score + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_analysis)
summary(rightwing1930_control_lm1)

rightwing1930_control_lm2=lm(rightwing1930 ~ BDEI_score2 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_analysis)
summary(rightwing1930_control_lm2)

rightwing1930_control_lm3=lm(rightwing1930 ~ BDEI_score3 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_analysis)
summary(rightwing1930_control_lm3)

rightwing1930_control_lm4=lm(rightwing1930 ~ BDEI_score4 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_analysis)
summary(rightwing1930_control_lm4)

rightwing1930_control_lm5=lm(rightwing1930 ~ BDEI_score5 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_analysis)
summary(rightwing1930_control_lm5)

### CLUSTERED STANDARD ERRORS

rightwing1930_lm1.vcovCL=cluster.vcov(rightwing1930_lm1, data_weimar_fin_analysis$wkr)
rightwing1930_lm1.clse=coeftest(rightwing1930_lm1, rightwing1930_lm1.vcovCL)
rightwing1930_lm1.clse

rightwing1930_lm2.vcovCL=cluster.vcov(rightwing1930_lm2, data_weimar_fin_analysis$wkr)
rightwing1930_lm2.clse=coeftest(rightwing1930_lm2, rightwing1930_lm2.vcovCL)
rightwing1930_lm2.clse

rightwing1930_lm3.vcovCL=cluster.vcov(rightwing1930_lm3, data_weimar_fin_analysis$wkr)
rightwing1930_lm3.clse=coeftest(rightwing1930_lm3, rightwing1930_lm3.vcovCL)
rightwing1930_lm3.clse

rightwing1930_lm4.vcovCL=cluster.vcov(rightwing1930_lm4, data_weimar_fin_analysis$wkr)
rightwing1930_lm4.clse=coeftest(rightwing1930_lm4, rightwing1930_lm4.vcovCL)
rightwing1930_lm4.clse

rightwing1930_lm5.vcovCL=cluster.vcov(rightwing1930_lm5, data_weimar_fin_analysis$wkr)
rightwing1930_lm5.clse=coeftest(rightwing1930_lm5, rightwing1930_lm5.vcovCL)
rightwing1930_lm5.clse

rightwing1930_control_lm1.vcovCL=cluster.vcov(rightwing1930_control_lm1, data_weimar_fin_analysis$wkr)
rightwing1930_control_lm1.clse=coeftest(rightwing1930_control_lm1, rightwing1930_control_lm1.vcovCL)
rightwing1930_control_lm1.clse

rightwing1930_control_lm2.vcovCL=cluster.vcov(rightwing1930_control_lm2, data_weimar_fin_analysis$wkr)
rightwing1930_control_lm2.clse=coeftest(rightwing1930_control_lm2, rightwing1930_control_lm2.vcovCL)
rightwing1930_control_lm2.clse

rightwing1930_control_lm3.vcovCL=cluster.vcov(rightwing1930_control_lm3, data_weimar_fin_analysis$wkr)
rightwing1930_control_lm3.clse=coeftest(rightwing1930_control_lm3, rightwing1930_control_lm3.vcovCL)
rightwing1930_control_lm3.clse

rightwing1930_control_lm4.vcovCL=cluster.vcov(rightwing1930_control_lm4, data_weimar_fin_analysis$wkr)
rightwing1930_control_lm4.clse=coeftest(rightwing1930_control_lm4, rightwing1930_control_lm4.vcovCL)
rightwing1930_control_lm4.clse

rightwing1930_control_lm5.vcovCL=cluster.vcov(rightwing1930_control_lm5, data_weimar_fin_analysis$wkr)
rightwing1930_control_lm5.clse=coeftest(rightwing1930_control_lm5, rightwing1930_control_lm5.vcovCL)
rightwing1930_control_lm5.clse



###########################
### NSDAP AND DNVP COMBINED VOTE SHARE (JUL. 1932)
### TABLE A25 IN THE SUPPLEMENTARY MATERIAL
###########################

rightwing1932_lm1=lm(rightwing1932j ~ BDEI_score, data=data_weimar_fin_1932_analysis)
summary(rightwing1932_lm1)

rightwing1932_lm2=lm(rightwing1932j ~ BDEI_score2, data=data_weimar_fin_1932_analysis)
summary(rightwing1932_lm2)

rightwing1932_lm3=lm(rightwing1932j ~ BDEI_score3, data=data_weimar_fin_1932_analysis)
summary(rightwing1932_lm3)

rightwing1932_lm4=lm(rightwing1932j ~ BDEI_score4, data=data_weimar_fin_1932_analysis)
summary(rightwing1932_lm4)

rightwing1932_lm5=lm(rightwing1932j ~ BDEI_score5, data=data_weimar_fin_1932_analysis)
summary(rightwing1932_lm5)

rightwing1932_control_lm1=lm(rightwing1932j ~ BDEI_score + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_1932_analysis)
summary(rightwing1932_control_lm1)

rightwing1932_control_lm2=lm(rightwing1932j ~ BDEI_score2 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_1932_analysis)
summary(rightwing1932_control_lm2)

rightwing1932_control_lm3=lm(rightwing1932j ~ BDEI_score3 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_1932_analysis)
summary(rightwing1932_control_lm3)

rightwing1932_control_lm4=lm(rightwing1932j ~ BDEI_score4 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_1932_analysis)
summary(rightwing1932_control_lm4)

rightwing1932_control_lm5=lm(rightwing1932j ~ BDEI_score5 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_weimar_fin_1932_analysis)
summary(rightwing1932_control_lm5)

### CLUSTERED STANDARD ERRORS

rightwing1932_lm1.vcovCL=cluster.vcov(rightwing1932_lm1, data_weimar_fin_1932_analysis$wkr)
rightwing1932_lm1.clse=coeftest(rightwing1932_lm1, rightwing1932_lm1.vcovCL)
rightwing1932_lm1.clse

rightwing1932_lm2.vcovCL=cluster.vcov(rightwing1932_lm2, data_weimar_fin_1932_analysis$wkr)
rightwing1932_lm2.clse=coeftest(rightwing1932_lm2, rightwing1932_lm2.vcovCL)
rightwing1932_lm2.clse

rightwing1932_lm3.vcovCL=cluster.vcov(rightwing1932_lm3, data_weimar_fin_1932_analysis$wkr)
rightwing1932_lm3.clse=coeftest(rightwing1932_lm3, rightwing1932_lm3.vcovCL)
rightwing1932_lm3.clse

rightwing1932_lm4.vcovCL=cluster.vcov(rightwing1932_lm4, data_weimar_fin_1932_analysis$wkr)
rightwing1932_lm4.clse=coeftest(rightwing1932_lm4, rightwing1932_lm4.vcovCL)
rightwing1932_lm4.clse

rightwing1932_lm5.vcovCL=cluster.vcov(rightwing1932_lm5, data_weimar_fin_1932_analysis$wkr)
rightwing1932_lm5.clse=coeftest(rightwing1932_lm5, rightwing1932_lm5.vcovCL)
rightwing1932_lm5.clse

rightwing1932_control_lm1.vcovCL=cluster.vcov(rightwing1932_control_lm1, data_weimar_fin_1932_analysis$wkr)
rightwing1932_control_lm1.clse=coeftest(rightwing1932_control_lm1, rightwing1932_control_lm1.vcovCL)
rightwing1932_control_lm1.clse

rightwing1932_control_lm2.vcovCL=cluster.vcov(rightwing1932_control_lm2, data_weimar_fin_1932_analysis$wkr)
rightwing1932_control_lm2.clse=coeftest(rightwing1932_control_lm2, rightwing1932_control_lm2.vcovCL)
rightwing1932_control_lm2.clse

rightwing1932_control_lm3.vcovCL=cluster.vcov(rightwing1932_control_lm3, data_weimar_fin_1932_analysis$wkr)
rightwing1932_control_lm3.clse=coeftest(rightwing1932_control_lm3, rightwing1932_control_lm3.vcovCL)
rightwing1932_control_lm3.clse

rightwing1932_control_lm4.vcovCL=cluster.vcov(rightwing1932_control_lm4, data_weimar_fin_1932_analysis$wkr)
rightwing1932_control_lm4.clse=coeftest(rightwing1932_control_lm4, rightwing1932_control_lm4.vcovCL)
rightwing1932_control_lm4.clse

rightwing1932_control_lm5.vcovCL=cluster.vcov(rightwing1932_control_lm5, data_weimar_fin_1932_analysis$wkr)
rightwing1932_control_lm5.clse=coeftest(rightwing1932_control_lm5, rightwing1932_control_lm5.vcovCL)
rightwing1932_control_lm5.clse


